%%%% Security Rate Maximization With AN Code
%%%% Probles Formulation
%%%%   max_(Q,m) log det(I + H1*Q*H1^H/Noise)- log_det(I + H2*Q*H2^H/(H2*A*H2^H+ Noise))
%%%%         s.t. tr(Q) <= SNR   Q>=0

clc;clear;
Nt = 6;
Nr = 1;
Ne = 4;
noisePower = 1;
SNRTotaldB = -15:5:30; %%%%最大发射功率单位dB
SNRTotal = 10.^(SNRTotaldB./10);  %%%% 最大发射功率约束，单位dB
 
channelNum = 10; %%% 生成10次信道
[hBsUserSet] = channelGenerate(Nr,Nt,channelNum);%生成100次信道HBsUe，H维度NrxNt
[hBsEveSet] = channelGenerate(Ne,Nt,channelNum);%生成100次信道HBsEve，H维度NexNt

rateSet = zeros(channelNum,length(SNRTotal));%生成最终速率数组
rateSetCmp = zeros(channelNum,length(SNRTotal));%生成最终速率数组




for indexh = 1:channelNum
    hBsUser = hBsUserSet(:,:,indexh); %具体单次信道hBsUser实现
    ANull =  null(hBsUser) ;          % 信道hBsUser的零空间求取 
    %人工噪声协方差矩阵应置于信道hBsUser的零空间，例如 hBsUser*A =0  / A = ANuLL*X
    % X为A的替代优化变量
    dimA = size(ANull, 2);
    hBsEve = hBsEveSet(:,:,indexh); %具体单次信道hBsEve实现
    objCVXCollect = [];
    objCmpUnit = [];
    for indexp = 1:length(SNRTotaldB)
        %Q维度为NtxNt，qCmpUnit是对比方案，代表在各个天线上分配相同发射功率
        QCmpUnit = diag(SNRTotal(indexp)/Nt*ones(1,Nt));
        Xlower = zeros(dimA,dimA);
        Alower = ANull * Xlower * ANull';
        %%%% 计算基于协方差矩阵QCmpUnit的安全可达造率
        rateCmpUnit = 1/log(2)*log_det(eye(Nr)+hBsUser*QCmpUnit*hBsUser')...
            - 1/log(2)*log_det(eye(Ne)+hBsEve*QCmpUnit*hBsEve'*   inv(hBsEve*Alower*hBsEve'+ noisePower*eye(Ne))  );

        objCmpUnit =  [objCmpUnit, max(rateCmpUnit,0)];
        Qlower = QCmpUnit;
        %%基于一阶泰勒展开的安全速率下界
        %log_det(noisePower*eye(Ne)+hBsEve*Q*hBsEve'+ hBsEve*A*hBsEve") <=..
        %log_det(sumLower)...
        %- real(trace(hBsEve"*inv(sumLower)* hBsEve*Qlower))...
        %- real(trace(ANull'*hBsEve'*inv(sumLower)* hBsEve * ANull*Xlower )).
        %+ real(trace(hBsEve"*inv(sumLower)* hBsEve*Q))...
        %+ real(trace(ANull'*hBsEve'*inv(sumLower)* hBsEve * ANull * X ))
        objBefore =10;
        objAfter = 10 * objBefore;
        count = 0;
        objAfterSet = [];
        while(abs((objAfter-objBefore)/objBefore)>= 5e-4 && count~=50)
            fprintf('第%g个信道实现，SNR为%g时的,第%g次内部迭代\n',indexh, SNRTotaldB(indexp), count);
            objBefore = objAfter;
            %%%%%%% CVX 使用框架
            cvx_begin quiet
            variable Q(Nt,Nt) hermitian semidefinite
            variable X(dimA, dimA) hermitian semidefinite
            maximize  log_det(eye(Nr)+hBsUser*Q*hBsUser') + log_det(  noisePower*eye(Ne) + hBSEve*ANull*X*ANull'*hBsEve')...
                      - real(trace(hBsEve'*inv(sumLower)*hBsEve*(Q + ANull*X*ANull')  ))
            subject to
            real(trace(Q + ANull * X * ANull')) <= SNRTotal(indexp);
            cvx_end
            %%%%%
            if(cvx_status(1)=='S')
                Qlower = Q;  % 更新Q，实现迭代的下界最大化
                Alower = ANull*X*ANull';
                objAfter = 1/log(2)*log(det(eye(Nr)+hBsUser*Q*hBsUser'))...
                         -1/log(2)*log(det(eye(Ne)+hBsEve*Q*hBsEve'* inv(hBsEve*Alower*hBsEve' + noisePower*eye(Ne)) ));  %CVX优化后的rate值
            end
            count = count + 1;
            objAfterSet=[objAfterSet, objAfter];
        end
        objCVX = max(objAfter,0);
        objCVXCollect = [objCVXCollect, objCVX];
    end
    rateSet(indexh, :) = objCVXCollect;
    rateSetCmp(indexh, :) = objCmpUnit;
    fprintf('The %gth channel realization is over \n', indexh)
end
